require(Rcpp)
output_version <- '2020-09-20'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,st_drop_geometry(ntl.zones))

cmr_prj <- read_sf(paste0(wkdir,'/local/gis_data/003_boundaries/cmr'),layer='localite_cmr')
cmr <- st_transform(cmr_prj,crs=4326)
ner <- read_sf(paste0(wkdir,'/local/gis_data/003_boundaries/ner'),layer='ner_political_cities')
nga <- read_sf(paste0(wkdir,'/local/gis_data/003_boundaries/nga'),layer='nga_political_cities')
tcd <- read_sf(paste0(wkdir,'/local/gis_data/003_boundaries/tcd'),layer='tcd_cap_admin')

cmr_join <- st_join(cmr,subset(ntl.zones,select=c(OBJECTID,GID_0)))
ner_join <- st_join(subset(ner,select=c(DESCRIPTIO)),subset(ntl.zones,select=c(OBJECTID,GID_0)))
nga_join <- st_join(nga,subset(ntl.zones,select=c(OBJECTID,GID_0)))
tcd_join <- st_join(subset(tcd,select=TYPE),subset(ntl.zones,select=c(OBJECTID,GID_0)))

##
cmr_join_out <- st_drop_geometry(cmr_join) %>% as.data.frame()
save.dta13(cmr_join_out, paste0(wkdir,"/proc_data/POLIcmr",output_version,".dta"))
ner_join_out <- st_drop_geometry(ner_join) %>% as.data.frame()
save.dta13(ner_join_out, paste0(wkdir,"/proc_data/POLIner",output_version,".dta"))
nga_join_out <- st_drop_geometry(nga_join) %>% as.data.frame()
save.dta13(nga_join_out, paste0(wkdir,"/proc_data/POLInga",output_version,".dta"))
tcd_join_out <- st_drop_geometry(tcd_join) %>% as.data.frame()
save.dta13(tcd_join_out, paste0(wkdir,"/proc_data/POLItcd",output_version,".dta"))

require(nngeo)

## CMR
# natcap
ntl.zones.pts["dist2natcap_cmr"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(cmr_join,Type=="Chef lieu national"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm1cap_cmr"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(cmr_join,Type=="Chef lieu regional"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm2cap_cmr"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(cmr_join,Type=="Chef lieu regional"), k = 1, returnDist = T)$dist) / 1000

## NER
# natcap
ntl.zones.pts["dist2natcap_ner"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(ner_join,DESCRIPTIO=="National Capital"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm1cap_ner"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(ner_join,DESCRIPTIO=="Admin1 Capital"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm2cap_ner"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(ner_join,DESCRIPTIO=="Admin2 Capital"), k = 1, returnDist = T)$dist) / 1000

## NGA
# natcap
ntl.zones.pts["dist2natcap_nga"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(nga_join,CLASS==1), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm1cap_nga"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(nga_join,CLASS==2), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm2cap_nga"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(nga_join,CLASS==3), k = 1, returnDist = T)$dist) / 1000

## TCD
# natcap
ntl.zones.pts["dist2natcap_tcd"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(tcd_join,TYPE=="Chef lieu national"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm1cap_tcd"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(tcd_join,TYPE=="Chef lieu regional"), k = 1, returnDist = T)$dist) / 1000
# regcap
ntl.zones.pts["dist2adm2cap_tcd"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(tcd_join,TYPE=="Chef lieu departemental"), k = 1, returnDist = T)$dist) / 1000

save.dta13(st_drop_geometry(ntl.zones.pts), paste0(wkdir,"/proc_data/POLIDIST",output_version,".dta"))
